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There are many instances in nature where geometry of the system that sustains chemical reactions 
is structured and the living cell is a typical example. There is a high degree of compartmentalization 
in the living cell where various compartments sustain different chemical reactions and transport 
reactants among themselves. In order to avoid storage facilities reactants are routed in orderly 
fashion between various places in the cell interior with large degree of synchronization. It is exactly 
such situation that is investigated here. The main motivation behind this study is to understand 
interplay between reactions and transport in a geometries that are not compact. Typical examples 
of compact geometries are box, sphere, etc. On the other hand, a network made of containers 
Ci, C2, • • • , Cjv and tubes is a typical example of space that is structured, and such non compact 
space is main focus on this study. The whole space is divided into a two regions. First, in containers 
particles react with rate A. Second, tubes connecting containers allow for exchange of chemicals 
with transport rate D. In such a way the two most important processes are isolated in the problem, 
reactions and transport. By varying topology of such network and details of chemical reactions it 
is possible to gain some understanding of interplay between chemical reactions and transport in 
structured spaces. It assumed that the number of reactants in the system is so small that kinetics 
is noise dominated. Two methods for solving corresponding master equation are discussed. The 
computer simulation is easy to implement, but leads to the results that are not that accurate. In 
here, a method is presented that can be used to calculate average, variance, and higher moments 
of the time reaction needs to finish. The method relies on a matrix representation of the master 
equation and is in principle exact. It works for an arbitrary reaction scheme and network topology. 
A number of different chemical reactions were studied and their performance compared in a various 
ways. Reactions are grouped into two ensembles, Reaction on a Fixed Geometry Ensemble (ROGE) 
and Geometry on a Fixed Reaction Ensemble (GORE). The ROGE and GORE are used to classify 
reactions in order to gain some understanding of which types of chemical reactions draw most benefit 
from the structured spaces. Most important findings are as follows, (i) There is a large number of 
reactions that run faster in a network like geometry. Such reactions contain antagonistic catalytic 
influences in the intermediate stages of a reaction scheme that is best dealt with in a network like 
structure, (ii) Antagonistic catalytic influences are hard to identify since they are strongly connected 
to the pattern of injected molecules (inject pattern) and depend on the choice of molecules that have 
to be synthesized at the end (task pattern), (iii) The reaction time depends strongly on the details 
of the inject and task patterns. 



I. INTRODUCTION 

The goal of the present study is to impact some 
progress in understanding interplay between reactions 
and transport in structured spaces. There are many ques- 
tions one could ask. This study focuses on one: How 
strongly the geometry of the system (e.g., topology and 
shape of the boundaries) influences the particular reac- 
tion schemes it harbors? This question can be asked 
in two ways. First, assuming that spatial structure is 
fixed, which reaction scheme would draw most benefit 
from it, for example, in terms of speed (smaller execu- 
tion time) or better timing (reduction of noise level)? 
Second, given the chemical reaction, what is the topol- 
ogy that is mostly suited for it? In order words, in which 
ways should one alter geometry of the system (without 
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altering reactants or influencing reaction mechanisms in 
any other way) in order to speed up chemical reaction 
encapsulated inside? Answers to these questions are rel- 
evant to a number of topics ranging from understanding 
of physio-chemical processes in the living cell, biological 
evolution, towards chemical engineering and biotechnol- 
ogy 

Taking the living cell as an example. The living cell 
exhibits large degree of spatial organization. Various 
compartments sustain different chemical reactions and 
transport reactants among themselves, with a high de- 
gree of spatial and temporal synchronization. Reac- 
tants and products of reaction are transported where 
they are needed, exactly at the right place and at the 
right time. There is no need of excessive storage of 
molecules. For reviews on these topic see 0, Q, and 

for other interesting 
work. The topics related to temporal and spatial syn- 
chronization of the cellular machinery are still hotly de- 
bated and present work might offer some understanding 
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of these phenomena. Furthermore, thinking from evolu- 
tionary perspective, one wonders what is that cells try to 
optimize by adopting highly structured geometry? 

Evolutionary process occurred in two stages: the chem- 
ical evolution |13|. which lead to emergence of organic 
molecules necessary for synthesis of proteins and any- 
thing alive, followed by biological evolution EiL a pro- 
cess which started after emergence of the first living cell. 
Fairly little has been done in understanding the role of 
geometry and spatial organization played in the evolu- 
tion process. There is already some pioneering work in 
this area 

El El El El E! Present study could impact 
some progress in this direction. 

Besides the topics discussed, the question of interplay 
between geometry and reaction scheme is relevant to a 
number of disciplines. The chemical engineering and 
biotechnology are typical examples. In general, there 
is a tendency to move away from bulk situation where 
volumes are large and reactants many, towards more ex- 
otic regime where reaction volumes are small and struc- 
tured and number of reactants relatively few. The inter- 
esting experimental work on these topics can be found 
in refs. [H El El E3 and El ^ The first set 
of references deals with reactions in liposome networks, 
and second set addresses the situation where a number 
of perfectly stirred reactors is connected by tubes. The 
exchange of chemicals through tubes is driven by pumps. 
The work done in here is relevant to both of these studies. 

The particular model proposed here is influenced by 
several lines of research. The setup of the reaction 
schemes is inspired by the work done on prebiotic evolu- 
tion, genome-based evolution and random reaction net- 
works liiiiii I^J^L and diffusion- 
controlled reactions in bulk pha se Ip , JM IH IM IsfllEcj 
and in restricted geometries [ill l42l liil l44| . The geomet- 
rical setup of the model is motivated b y th e experimental 
studies given in [H El El El El |H El and theo- 
retical stud y o f reaction-diffusion neuron and enzymatic 
neuron [H IH El El El El E2 

The situation considered here is radically different from 
the bulk studies. The goal is to mimic particular as- 
pect of the cell environment where reactions happen in 
the specific regions of space and reactants move among 
these regions. The simplest way of achieving such a setup 
is to consider chemical reactions in containers that are 
connected by tubes. Figure ^shows one possible exam- 
ple. Reactants can move from container to another along 
tubes. Thus, the compartmcntalization is build into the 
model in the simplest possible way. Furthermore, it is as- 
sumed that containers are small in size so that the num- 
ber of reactants in each container is also small. Such 
setup inevitably leads to fluctuations, and chemical ki- 
netics becomes noisy, as is usual the case in the living 
cell. 

For simplicity reasons, the details, both of the chemical 
reactions and transport, are neglected to a large extent. 



It is assumed that reactants are point like objects with no 
structure. What is studied here can be classified as arti- 
ficial chemistry . The two most important time scales 
are traced in the model. The reaction rate A describes 
how fast molecules react in containers. The transport 
rate D governs exchange of chemicals among contain- 
ers. It is assumed that diffusion is the dominant mech- 
anism of particle transport and this is reasonably good 
approximation for the living cell. Q, Naturally, there 
are exceptions to this rule, e.g. transport of molecules 
by kinesin along microtubule in the cell, but such cases 
are not going to be considered here. For discussion of 
reaction-diffusion systems in biology please see ref. j5^| ■ 

How to compare chemical reactions in a compact and 
structure geometries depicted in Fig. [panels (a) and 
(b)]? Naturally, the most obvious way would be to run 
chemical dynamics directly on structures depicted in pan- 
els (a) and (b) and try to trace differences. However, it 
will be shown that this can also be done by studying ki- 
netics in network like geometry depicted in panel (c), but 
such strategy has to be implemented carefully. For exam- 
ple, it is easy to see that panel (c) in Fig.^is a very rough 
representation of structured geometry in panel (b). How- 
ever, there is no apparent similarity between structures 
shown in panels (a) and (c). Nevertheless, the network 
structure in (c) contains (a) and (b) as special cases. For 
example, when D ~ A the reaction dynamics occurring in 
the network depicted in panel (c) should exhibit roughly 
the same features as the reaction dynamics running in 
the structured geometry depicted in panel (b). Further- 
more, when transport rate dominates all other processes 
in the system, D ^> A, same should hold for (a) and (c). 
Also, in such a case, D should be much larger than any 
other reaction rate that might enter due to the effects 
of catalysis. Thus, in order to gain some understand- 
ing of differences between (a) and (b) one simply has to 
study variations of a reaction dynamics in the network 
like structure (c) as transport rate changes from D ~ A 
towards D 3> A. 

The type of the analysis used is largely static. To 
understand in which ways particular shape of the reac- 
tion volume influences chemical kinetics (and vise versa) 
a large number of chemical reactions will be analyzed 
and their performance compared. The set of reactions 
and geometries chosen form an ensemble A — {o w ,w = 
1, . . . , E} where elements in ensemble are pairs o u = 
(reaction, geometry). The element o w will be referred to 
as an organism. 

Two ensembles are distinguished. In the first case, the 
reaction scheme is kept fixed, while the geometry of the 
network is allowed to change. This type of ensemble will 
be referred to as the geometry on the (fixed) reaction en- 
semble (GORE). The geometries are different and they 
are sampled by changing size of the containers and length 
of the tubes. In the second case the geometry of the net- 
work is kept fixed, while the reaction scheme is subject to 
a change. This type of the ensemble will be referred to as 
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the reaction on the (fixed) geometry ensemble (ROGE). 
In this work most of the attention is on the ROGE en- 
semble. 

What is the measure of the good performance for a 
chemical reaction? In here, focus will be on the time re- 
lated issues such as the length of catalytic cycles. The 
importance of the timing in the intra cellular process has 
been discussed by Volkenstein 'The most impor- 

tant role in ontogeny is played by time factors, i.e. the 
time relations in the synthesis of various proteins and 
the formation of various cells and tissues. Even small 
changes in the timetable lead to considerable morpholog- 
ical alternations...'. The setup presented here allows for 
consideration of other performance criteria but these are 
omitted due to the simplicity reasons. 

The time needed for a reaction to finish, T, is a 
stochastic variable and for a given reaction scheme and 
topology it will have well defined distribution function 
<!>(£): the probability that reaction finished within the 
time interval [t, t + dt] is given by P(t < T < t + dt) = 
<fr(t)dt. All information about reaction timing is hidden 
in the distribution function <fr(t). In practise, it is very 
hard to find $(£) for general reaction scheme and geome- 
try. Also, amount of information contained in is too 
large. It is more fruitful to characterize Q(t) in terms of 
few variables and in this work we use two: the average 
time needed for reaction to complete r = J °° t$(t)dt and 
the variance a where a 1 — J* °°(i — r) 2 ^(t)dt. In such a 
way it is possible to asses how fast reaction happens and 
how noisy it is. 

Once performance criteria are quantified one faces op- 
timization problem of finding the best performer in a 
given ensemble A. In principle, the best performer will 
be the organism that exhibits smallest value for r. Two 
strategies are used to find such organism. First, one 
simply generates full ensemble, measures performance of 
each element, and sorts performance at the end. How- 
ever, when number of elements in such ensemble is large 
another optimization method has to be used. For ex- 
ample, when there is a need for that , on e can use a 
method described in refs. [H |U |H M, E3, E3 
This optimization method mimics process of evolution 
and belongs to the class of genetic optimization algo- 
rithms. The terms evolution and optimization will be 
used interchangeably through out the text. 

The paper is organized as follows. In section [D] the 
reaction-diffusion model is defined in detail, and in sec- 
tion II I II overview of the methods for solving the master 
equation is given. The method of doing a computer sim- 
ulation is reviewed first. The method of finding moments 
of the time length of the catalytic cycles is discussed in 
detail in the next section llVl Few simple reactions are 
discussed in section where it is shown how to compare 
vastly different reaction schemes. Section IVII contains 
analysis of ROGE ensemble made from two container 
network, two particle types, all possible reactions, and 



full set of inject-task patterns. Conclusions and outlook 
are given in section IV 111 



II. THE MODEL 

The reaction-diffusion model is defined as follows. 
Consider the set of containers Ci connected in a partic- 
ular way by tubes with lengths U^; i,j = 1, . . . , N. It is 
possible that some containers are not connected. Exam- 
ple of such structure is given in Fig. 1. The containers 
harbor molecules X a ; a = 1,2, ... ,M. Molecules are 
allowed to react only when in the same container (pro- 
vided there is a reaction they participate in). Molecule 
X a moves from the container Ci to the container Cj with 
a rate D"" . Please note that the expression for the trans- 
port rate is diagonal in the particle index a (there is no 
change of particles while being transported). This as- 
sumption could be easily relaxed. For simplicity reasons 
it is assumed that Df" = f(kj) if the link between con- 
tainer i and j exists, otherwise Dff = 0. Links are non- 
directional and transport rates are same for all particles, 
i.e. D™™ = Dij = Dji. Please note that the regular 
lattice is a special case of this very general model. Struc- 
tures like these can be found in the interior of the living 
cells. Mapping is not exact, but there is a rough similar- 
ity. 

By assumption, the reactions only occur in the con- 
tainers and not in the tubes, and reactants mix well in 
the containers. Tubes serve only as connectors of regions 
where reactions happen. With assumptions at hand, to 
describe system at any time instant, it is sufficient to 
track number of particles in each container. This greatly 
simplifies calculations. Conformation of the system c is 
specified as a occupancy of the containers, 

c= (m,. ..,ni,..., tijv) (1) 

where vectors rii i — 1, N describe particle content of 
each container, 

rii = (ni,i,...,n a}i ,...,nM,i) (2) 

with n a ^ = 0, 1, 2, oo for a = 1, M and i = 1, N. 
System makes random transitions between various con- 
figurations. A configuration of the system, c, changes 
when either reaction or transport occur. 

For computational convenience, the very simple type 
of a reaction scheme will be considered, 

X a -> T Xp (3) 

where X 1 is understood to influence conversion of X a to 
Xp either as catalyst (+) or suppressor (— ). All reac- 
tions that are allowed are assumed to have same reaction 
rate A. The reaction graph is specified by the reactivity 
matrix A. If reaction X a — > Xp is allowed A Qj/ 3 = 1 and 
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equals zero otherwise. The reaction graph is directed and 
matrix A is not necessarily symmetric. 

The effect of catalyst or suppressor are modeled as fol- 
lows. The bare reaction rate A may be modified by the 
presence of other reactants in the container C,;. It is as- 
sumed that each reaction can have at most one catalyst 
or suppressor. On the other hand, it is possible that one 
particle is catalyst (or suppressor) for more than one re- 
action. For programming purposes it is sufficient to use 
array K ai p = ±7 if X 1 is catalyst (+) or suppressor (— ) 
for X a — > reaction. K a ,p — indicates that reac- 
tion does not have catalyst nor suppressor. Also, it is 
assumed that the effects of the catalysis are strongly en- 
hanced. If there is some catalysis going on, it completely 
dominates reaction scheme. Using these assumptions the 
reaction rate for the reaction X a — > Xp in the container 
Ci is given by 

f 1 k = 

Kx,p( ni ) = XAa -P { Z K>0, 71k - < > SK ' a ( 4 ) 
[ i K < 0, n\ K \ ti > 5\ K \ <a 

S a> /3 is a Kronecker delta-function, £ > 1 denotes catalysis 
enhancement factor, and k = K a p. Equation (0J is self 
explanatory for the exception of possibly one term that 
we proceed to discuss. The condition n Kt i > 5 K , a ensures 
that there is no self catalysis for the reaction of the type 

X a ^ X/3 (5) 

when there is only one X a present in the container. One 
needs at least two X a in order to fell catalytic influence 
of X a on the reaction given in JSJ . 

The particular type of the reaction schemes considered 
here (0) is inspired by chemical processes in the living 
cell. It is often that by action of enzymes certain input 
set of chemicals is converted into output molecules by se- 
ries of intermediate reactions. It is exactly this aspect of 
cellular machinery which this scheme tries to capture in 
the most simple way. The most general reaction graph 
one can consider is shown in Fig. [21 The graph indicates 
that cell machinery converts molecules Xi,...,X u into 
molecules X[ , X' 2 , . . . , X^ . The shaded area in the mid- 
dle denotes intermediate reaction steps that involve ar- 
bitrary set of reactions between molecules already shown 
in the graph. An additional particle types may appear 
in the shaded region. It is possible to define the speed 
of a reaction as the time needed for the predetermined 
set of output molecules X'^X'2, ...,X' n to appear for the 
first time, provided only the input molecules were present 
initially in the system. How this idea is implemented is 
shown bellow. 

In addition to reaction scheme given in Eq. (|3J) two 
quantities are specified. First, at t = a certain number 
of particles is injected into the system in various contain- 
ers. The list of the particles injected is specified by the 
vector 



In the course of time the injected set of particles will be 
transformed into something else. Second, the vector 

7T = (n* 1A , ...,<,,... ,n* MN ) (7) 

specifies tasks that have to be achieved. For example, 
n* j > indicates that goal is to synthesize n* a i molecules 
of the type X a in the container Ci . On the other hand, 
n a i — indicates that the number of particles X a in 
the container Ci is not traced. The maximum number 
of tasks is given by N x M, though not all tasks need 
to be monitored. Also, the cases where task is achieved 
trivially at t = are forbidden and additional restriction 
is set upon components of 1 and 7r: If task is traced one 
has n* ai > 0. In such case, the condition n* ai > nf\ 
must hold. 

To make a notation easier, it is convenient to eliminate 
the elements from l and ir that are zero and write 

t — 0*1 , ^2, ■ ■ ■ , O, n — C 71 "!) n 2' ■ ■ ■ > n v) (^) 

where /., and 7r, only contain list of molecules injected, 
and tasks that are actually monitored. It is clear that 
uj,ri < M x N. 

Every time a certain task is accomplished the time 
when this happens is stored. These times are arranged 
in the vector 

T=(T 1 ,T 2 ,...,T V ) (9) 

and there is a one to one correspondence between the el- 
ements of 7r and T. Once the task is achieved molecules 
that were used to accomplish it are removed from the 
system. This consideration is motivated by the character 
of the real processes in the living cell. If certain number 
of molecules are needed at specific place, once arriving 
there, these molecules will be consumed by other bio- 
chemical processes in the living cell. 

The vector T is a stochastic variable and can 
be described in terms of the distribution function 
&(ti, t-2, . . . , t„; i, 7r). In practise, it is very hard to ob- 
tain full distribution function and it is more convenient 
to use first two moments, the average 

t = (n, -,T k , ...,t„) (10) 

and the variance 

f = (oi, 1) •••) Oa.i) •••> O-JW.JV) (11) 

In addition, one could also include moments of the type 
(rjTj) i ^ j but these will be not considered. 

The quadruple consisting of particular reaction scheme 
(A, £, A and K), network geometry (l iy j i,j = l,..., N), 
inject pattern t, and list of tasks monitored 7r, will be 
referred to as an organism. The organism can be seen 
as the entity that has to transform a certain number of 
chemicals into a set of molecules that have to be synthe- 
sized at certain places, utilizing available reaction scheme 
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and geometry. In the following sections various organ- 
isms will be classified according to the criteria how fast 
they achieve certain list of tasks (see. Ea. 110(1 . In prin- 
ciple, it is possible to make classification with regard to 
how noisy performance is (see Ea. Illj l. but this type of 
classification is left for future studies. 



III. THE DYNAMICS 

The dynamics of the system defined in the previous 
section is stochastic and from the rules discussed one can 
derive a master equation that describes the time evolu- 
tion of the occupation probabilities p(c, t), 

P(C, t) = Rc,a'P(c', t)-J2 R C ' CP {C 7 1) (12) 
c' c' 

where here and in the following dot over symbol denotes 
time derivative. The reaction rates R c 'c describing tran- 
sition c — > c' can be easily calculated from the definition 
of the model. In general, it is very hard to solve the 
equation i|12|) . In here, two strategies are use to solve it. 
First strategy is based on a simulation method. This is a 
straight forward approach. In the second instance a set 
of equations is derived that specifies first, second, and 
possibly higher moments of the $(ti, . . . , t n ; i, it). Both 
strategies are implemented into a computer program. 

Computer Simulation: The master equation (I12|) is 
solved by using minimal process algorithm suggested by 
Gillespie [5J, |55| . Given that the system is in a certain 
configuration, one can calculate distribution of waiting 
time for the next process to happen. Time is updated 
by amount of waiting At. Process is chosen randomly 
according to weight given by corresponding rates for each 
process. For problem at hand a linear selection algorithm 
is used. 

Given that at the time t the system was in the confor- 
mation c specified in ft), following processes can happen. 
Transport of particle X a from Ci to Cj occurs with rates 

R? tj (c) = D? tj n a ,i , a = l,2,...,M, i,j = l,...,N 

(13) 

or reactions within containers with rates 

KA C ) = (14) 
One also needs the total reaction rate, 

N M N M 

Q( c )= E E^( c )+E E R «A c ) ( 15 ) 

i.j — 1 a—1 i—1 a,/3— 1 

which is used to calculate the time update At = — ln(l — 
r)/Q where r is a random number drawn uniformly from 
the interval [0, 1). The probabilities for process to happen 
are given by 

RfJc) 

/<:\ —fz-- a = l,...,M, i,j = l,....,N (16) 



Pa,/s = - 5 g iji , i=l,...,N, a,p = l,....,M(lT) 

A process is chosen using the linear selection algorithm. 
First, a new random number r' is drawn. After that, the 
cumulant probabilities are calculated through loop over 
all processes. The loop is stopped when the cumulant 
probability exceeds r' and the last process for which this 
happened is executed. 

Moment method: This method relies on the matrix 
representation of the master equation (|12fl and can not 
be used to treat cases with large number of containers and 
particle types. However, the method is exact and should 
be used when there are enough computation resources. 
The method is developed in the following section. 



IV. THE CALCULATION OF THE AVERAGE 
AND STANDARD DEVIATION 

It will be shown how to calculate moments of the in- 
dividual components of r, 7^ = ((ifc) p ), k — 1, . . . , 77 
and p = 0,1,2,. ..,00. The more complicated mo- 
ments, e.g. 7^'^ = {(Tk) p (n) q ) with k, I = 1, . . . ,rj and 
p,q = 0, 1, 2, ... , 00, could be also evaluated using tech- 
nique presented in this section, but such moments will 
not be considered. The focus in on the moments that are 
obtained through 

poo 

j { k p) (L,TT) = / t p T k (t;L,n) dt (18) 
Jo 

where Tk(t; t, 7r) denotes integrated distribution function 
for task tt/- given by 

r k (t;t,TT)= / ^(tx,...,t n ;L,w) Y[ dt m (19) 

Please note that accomplishment of each individual task 
is influenced by presence of others since the particles can 
vanish upon accomplishment of various tasks. This is the 
reason why integrated distribution function Tk{t; l,tt) 
contains both index of the task that statistics is sought 
for (fc) and full list of tasks being monitored (ir). 

It is possible to find closed expression for Laplace 
transform of Tk(t;L,ir). The Laplace transform 
of arbitrary function F(t) is defined as F(s) = 
J exp(—st)F(t)dt. Tk(s; l, it) is given by 

r fc (s; t, 7r) = 2J w(c, ir k )g(s; 1, c, n) 

+ E E w ( c ' 7rm ) 5 ( s;t ' c ' 7r ) 

771=1,77 c ¥^ L 

xT k (s; c/7r m , 7r/7r m ) (20) 
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where notation used is as follows. The w(c,TTk) equals 
one if task can be accomplished once system arrives 
in the state c, and equals zero otherwise. In the follow- 
ing, by definition, a state for which one of the w(c,TTk) 
with k = 1, ... ,77 differs from zero will be referred to 
as a window state. Through window states tasks can be 
accomplished. c/n m denotes the state immediately af- 
ter the particles have been taken away once the task 7r m 
was accomplished. Likewise, the symbol 7r/7r m denotes 
a list of tasks being monitored with task 7r m omitted 
(7r 1 ,...,7r m _i,7r m+ i,...,7r t? ). g(s; t,c,7r) is a distribu- 
tion function for the first passage time into the state c 
given that the dynamics started from the state l. Fig. |3| 
is a schematic presentation of the Eq. Q20(l . Please note 
that Eq. (|20|) defines Tk(s;L,Tz) recursively. The num- 
ber of tasks on the right hand side of Eq. I|2U|I is by one 
smaller than the same number on the left hand side of 
the equation. When list of a tasks is empty tt — ttq and 
7Tq = ( ). Condition Tk(s; /,, 7T ) = stops recursion. 

The Laplace transform of the first arrival time distri- 
bution function g(t; l, c, tt) is calculated as follows. Given 
the particular reaction scheme it is possible to construct 
master equation l(T2"|) that governs the time dependence of 
the occupation probabilities of each state p(c, t) , where c 
has to be accessible from the initial state l. When calcu- 
lating matrix of the transition rates R it is assumed that 
all window-states can not be left once they are arrived 
into. Window states are perfectly absorbing. Fig. is a 
graphic representation of this fact. Once p(c,t) is found 
from (|12|) the first passage time distribution function is 
given by g(t, t, c, tt) = p(c, t). 

It is useful to arrange both p(c, t) and <?(s; t, c, 7r) into 
a vectors p(t) and g(t) where notation was simplified a 
bit since we assume that 1 and tt are known and fixed. It 
is useful to rewrite master equation (|12|l in a matrix form 
as p(t) = Rp(t). This equation is solved using Laplace 
transform, with initial condition p(c, 0) = 8 C L (5 denotes 
Kronecker delta symbol): sp(s) — Po = Rp(s). Also, in 
the Laplace transform space one has sp(s) — po = g(s), 
which directly leads to the equation for the first arrival 
time distribution function: 



sg(s) = R[g(s)+p Q ] 



(21) 



In principle, the equation above could be solved as g(s) = 
(s—R)~ 1 pq. However, matrix R has zero eigenvalues and 
the value of g(s) in the limit s — > is ill defined. To avoid 
such problems Eq. I|21|l has to be solved in a special way. 

It is useful to separate configuration space into three 
groups, as shown in Fig. [3] First group, labeled S n , con- 
tains states that are non-window or the normal-states. 
Second group contains states labeled by St that we re- 
fer to as the trap states. The third group consists of the 
window states solely, labeled by S w . The existence of the 
trap states is problem dependent. Once system arrives 
into these states there is no exit from this space, though 
such states are not window states. This simply means 
that it is possible that set of tasks is never accomplished. 



Using the partition of states shown in Fig. [3] leads to 
the following set of equations 

sg n {s) = Rnn[9n(s) +Pn,o] (22) 

sg t (s) = R tn [g n (s) +p n ,o] + Rttgt(s) (23) 

sg w (s) = Rwn[g n (s) +Pn,o] (24) 

Please note that pt,o an d Pw,o are zero since, initially, the 
system is in the state t and such state does not have any 
components in the St and S w spaces. Given that there 
are no transition from trap states into normal states or 
window states blocks R nt and R w t are missing in the 
equations above. Likewise, blocks R ww , R n w and Rt w 
and are zero since there are no transitions among window 
states, nor transitions from them. 

The solution of the equations (|22|l - 124|) can be found in 
a straight forward manner. Equation (|22fl can be solved 
first, leading to g n (s) = (s — Rnn) Pn,o an d inserting 
this expression into the Eq. (|24ll gives 

g(s, l, c,7r) = [R w „(s - R nn )~ 1 Pnfl] c , ceS w (25) 

Please note that Eq. i|25|) is well defined for all values 
of s. In particular, in the limit s — > even for a matrix R 
that has zero eigenvalues. It is intuitively clear that, con- 
trary to R, matrix R nn does not have zero eigenvalues: 
as time goes on, all probability accumulates in S w and 
St spaces (see Fig. |3J). The only difficulty with Eq. (f25|) 
is partitioning of the full configuration space into S n , S w 
and St- The algorithm for carrying out such partitioning 
is not presented here in order to save the space. 

Finally, once g(s, t, c, 7r) is found one can proceed with 
the calculation of the moments 7^(1, 7r). These can be 
obtained by taking derivatives of Eq. (|20|l with regard to 
s and setting s = at the end: it can be seen easily from 
Eq. O that 



7fc(i,w) 



(-)* Hm c?jr fc ( S , t ,7r) 

s— »0 



(26) 



where d s denotes derivative over s. Using Eqs. (12(i[) and 
PO) leads to 



(p) 



(t, 7T) = ^ W (°, 7Tfc)ff (P) (<-, C, TT) 



X 7fc P_9) ( c / 7r m I 7r / 7r m) 



(27) 



where by definition g( p ' (t, c, 7r) = 
(— ) p lim s ^ d^g(s, t, c, 7r), which after using Eq. lt25|) 
leads to 



<?W( t ,C,7T) 



(-) 



R n 
±tj wn ±tj nn fn.O 



(28) 



Equations (|27|I - H28|I are central result of this sec- 
tion. They determine all moments. For example, once 
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7^ (i, 7r) p = 0, 1, 2 are found the average and variance 
of T are given by 



7>,(t,7r) = 








(29) 










cr fc (t,7r) 2 = 






2 


(30) 


where /c = 1, . . . , r\. 


One has to divide by 7f°' ) (t 


7r) in 



the equations above in order to ensure that in the case 
when there is a possibility that some tasks are not ac- 
complished statistics is done only for instances where 
task was achieved. The percentage of cases when this 
happened is given by 7? (i., 7r). The numerical imple- 
mentation of Eqs. (|27(l and Ij28(l is a straight forward and 
gives the exact values for r and er. 



V. ROGE ENSEMBLE: DESCRIBING 
ORGANISM PERFORMANCE IN TERMS OF 
THE SINGLE VARIABLE v 

The calculation of t and er for an arbitrary organism 
was discussed in the previous section. In this section the 
GORE ensemble will be studied. The case N = 1 where 
there is only one container is not interesting since such 
space is not structured. The first non-trivial example of 
a structured space is the case of a two-container network 
with N — 2. For simplicity reasons, a situation will be 
considered where there are only two particle types A and 
B corresponding to M = 2. With only one particle type 
one can only focus on transport issues. To see coupling 
between reactions and transport one needs at least two 
particles types. Please note that the analysis done in 
this section is quite generic, though it is carried out on 
a rather simple case of N = 2 and M = 2. The analysis 
could be easily repeated for an arbitrary values of N and 
M. However, the computational cost scales with TV and 
M, and there is clearly an upper limit to which cases one 
can study. It will be shown that the relatively simple two- 
container network and reactions with two particle types 
provide interesting insight into the problem. 

The structure of the organisms in the ROGE ensemble 
is defined as follows. For each organism in the ensemble 
an unique choice is made for (i) total number of A and B 
molecules, (ii) reactivity matrix A, (iii) catalytic activ- 
ity matrix K , (iv) inject pattern 1 and (v) list of tasks 
monitored tt. For all organisms geometry is kept fixed 
(e.g. size of containers and length of the tube connect- 
ing them). Please note that there are two symmetries in 
the problem, the one that originates from relabeling of 
particles, and another one that has to do with relabeling 
of containers. A special care is taken to eliminate these 
symmetries in the ensemble. The goal is to unearth best 
reaction scheme (organism) that draws most benefit from 



the structured space that has the form of a two-container 
network. 

We start with the situation where the total number of 
particles N* = Ua,i + ha,2 + n B,i + nB,2 in the system 
equals one. Also, we consider only organisms with reac- 
tions A — > B and B — ► A, both with rate A. One can see 
that, for a given reaction scheme, and with a constraint 
N* = 1 , only three choices for l and 7r are possible, lead- 
ing to the three organisms o\ , 02 and 03 that are listed in 
table [I] Please note that all other N* — 1 cases can be 
related to these three through relabeling of particles and 
containers, or by considering different reaction schemes. 

All organisms are such that one A particle is injected 
in the container C\ and only one task is monitored. In 
the case of organism o\ the goal is to synthesize one A 
molecule in the container C2, in the case of 02 one B 
molecule should be synthesized in C\, while in the case 
of o 3 one B molecule should be created in Ci. Reaction 
process is stochastic. Average and standard deviation of 
time for task accomplishment are given under columns 
labeled r and a in table [I] The dependence of r and a 
on (the inverse of) the transport rate D is depicted in 
Fig.H 

Clearly, all three organisms achieve their tasks faster 
in the compact geometries. This can be seen from Fig.0] 
since r gets smaller when D _1 ,Z — » 0. (Here, we used 
the fact that D^ 1 ~ I 2 where I is the length of the tube 
connecting containers). Same holds for curves depicting 
a. Amount of noise decreases when geometry is compact. 
For example, organism o\ functions by transporting one 
A particle from container C± to container C%. It is intu- 
itively clear that this happens faster when containers are 
close. One can analyze 03 in the similar way. However, 

02 is somewhat different. In the case of 02 the goal is to 
synthesize one B molecule in the same container where A 
was injected. Figure 0] illustrates the fact that when an- 
other open volume is present molecule can wonder away 
into additional volume and this process delays synthesis 
of B. One can also see that 02 is least sensitive to increase 
(decrease) in I (D): when l.D^ 1 — > 00 the r for o\ and 

03 increases, while for 02 r saturates to constant (though 
level of noise, a, increases). 

The question is whether it is possible to present infor- 
mation conveyed from Fig. 0] and table Q] in a more com- 
pact way? There are couple of reasons for this. First, 
in a case of the more complicated geometry, generating 
graph such as shown in Fig.0]is time consuming. Second, 
it would be desirable to have automatic procedure for as- 
sessing most important properties of such graph; whether 
given organism performs best in compact or structured 
(network like) geometry. 

To achieve this goal and compactify information in 
Fig. 0] and table |I] it is useful to consider following quan- 



8 



tity, 



(31) 



where |.| denotes Euclidean norm of the vector, t„ and 
To are given by r calculated for network like and compact 
geometries with jump rates D n and Dq such that D n ~ 
A and Dq 3> D n ,\,£\. Using similar reasoning the v 
can be defined for generic network having more than two 
containers. In such a way one can compare extended and 
compact geometries of given fixed network structure and 
express comparison through one variable v. 

In the following the v will be refereed to as speed of a 
reaction, v < 1 is indication that organisms accomplishes 
tasks faster in a network like geometry, while v > 1 shows 
that organism draws most benefit from a compact geom- 
etry. Please note, all organisms considered in table|l]and 
Fig. U| have v > 1. 

In principle, one could define quantity similar to v and 
use a instead of r. For example, one could use 



(32) 



where meanings of cr n and ctq are similar to the ones of 
t„ and t q . [i could be used to classify organisms in the 
ensemble according to the amount of noise in compact 
and extended (network-like) geometries. Also, it is pos- 
sible to use p = y/ v 2 + fi 2 as simultaneous measure of 
speed and noise. But such quantities will not be studied 
at the moment. From now on we focus on v solely. 



VI. ROGE ENSEMBLE: CLASSIFICATION OF 
THE REACTION SCHEMES USING THE SPEED 
OF REACTION v 



Figure [3] depicts results of the classification of large 
number of organisms in five ROGE ensembles where v 
(defined in previous section) is used as a measure of 
the performance. Five ROGE ensembles are constructed 
with the increasing upper limit for the total number of 
particles in the system N* from 1 to 5. 

For N* — 1 there are 95 unique organisms. Clearly 
this number is overestimate. The number of unique or- 
ganisms is obtained after eliminating symmetries related 
to relabeling of particles and containers removed. How- 
ever, this number is still too large since catalytic influ- 
ences are assumed to play the role, though there can not 
be any catalytic influence when there is constantly one 
molecule in the system. The number of unique organ- 
isms for other cases discussed later with N* — 2, 3, 4, 5 is 
correct. There are no organisms in the N* = 1 class that 
benefit from the structured geometry since the histogram 
in the panel (b) is empty. All organisms have v > 1 as 
can be seen from panel (a). 



For N* = 2 there are 730 organisms that are unique. 
Only after more than one molecule appears in the system, 
catalytic influences start to play the role, and organisms 
that benefit from the network like structure appear [see 
panels (c) and (d) in Fig. [5]. The best performer in this 
class is given in Table ITT1 denoted by Oi.2, together with 
couple of a second best performers denoted by oa i- The 
two organisms labeled 0^2 draw most benefit from the 
network like structure and degeneracy in v comes from 
the fact that same enhancement and reduction factor is 
used for positive and negative catalytic influence respec- 
tively. 

One can understand intuitively why Oj.2 runs faster 
in the network like geometry. We focus on the particular 
case of Oi ; 2 where goal is to synthesize two B molecules in 
the container C\ . Due to the initial presence of molecule 
A in the container C\ , it is very likely that B molecule will 
be converted into A molecule, when A and B meet in the 
same container. Once there are two A molecules in the 
system the trouble starts. Even if the reaction A — > B 
happens, and chance for this is really small due to the 
negative catalytic influence of A on such reaction, B will 
be converted back to A immediately (due to the positive 
catalytic influence of another A molecule on the B — > A 
reaction). In principle, conversion of AB into 2B has no 
chance occurring in a reasonable time when there is only 
one container. The antagonistic catalytic influences just 
discussed are best handled in a network like geometry. In 
such a case the synthesis of the molecules can be done in a 
separate containers. With two containers there is always 
a chance that antagonistic influences will be reduced. For 
example, in the case there are two A molecules in the 
system, the processes A — > B can happen fast given that 
damage inflicting A molecule is in another container. 

Also, it is naive to think that organisms containing re- 
actions with solely negative catalytic influence perform 
best in a network like structure. This is clearly not the 
case. The winning organisms 0^2 are constructed from 

one reaction with positive B — > A and another reaction 

— A 

A — > B with negative catalytic influence. Furthermore, 
Fig. O shows plenty of other cases. Actually, the com- 
pletely opposite is possible. There are many organisms 
in the histogram plot with v > 1 that contain at least 
one reaction with negative catalytic influence and these 
organisms perform best in a compact geometry. 

With TV* = 3 there are 3025 organisms that are unique. 
The best performer in this class is given in Table UTI de- 
noted by Oi.3. A couple of the second best performers 
are also shown and labeled 0^3. When more than two 
molecules appear in the system, completely new organ- 
ism appears as winner. There is a sharp transition in 
character of winners. Both the reaction type and inject 
and task patterns are different in 0^3 and 0^2- The reac- 
tion of organisms 03^ is such that it is possible that tasks 
are not accomplished. For example, all A molecules can 
be converted to B before task [2A]-[ ] in C\ is accom- 
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plished. Once this happens synthesis of A's in C\ will 
never occur since back reaction B — > A is absent. 

Very likely, from biological point of view there is no 
advantage basing survival on a reaction that sometimes 
fails, i.e. using reaction similar to the one contained in 
0^3. Nevertheless, there might be constraints that could 
enforce presence of such reaction. For example, the num- 
ber of molecules in Nature is large but limited. Living 
organisms have to use what is available. In the lack of al- 
ternatives it might be necessary for the intra cellular ma- 
chinery to use reaction that sometimes fails. But leaving 
such discussions aside, the only point emphasized here is 
that organism 0^3 draws most benefit from the network 
like structure. 

Interestingly, this type of organism remains winner in 
the classes N* = 4 and N* = 5 with v = 0.0248 and 
v = 0.0192, see table|H] organisms 0^4 and 0^5. However, 
please note that details of the inject pattern of 0^4 and 
0^5 are somewhat different from 0^3. Judging solely from 
Oi t 3 and Oj,4 one would guess that inject patter for best 
performer in class N* = 5 should be t =[A] — [4A], in class 
N* = 6 t =[A]— [5A], etc. However, this is not the case. 
The best performer in class N* — 5 has inject pattern 
equal to l =[2A] — [3A] with v = 0.0192, while organism 
with l =[A] — [4A], denoted by 0^2, has somewhat larger 
v = 0.0219. 

Thus, examples above show how difficult it is to have 
any intuition about structure of best performers. There 
are also other ways to see this. For example, the last 
row in table [n] contains organisms with slightly modified 
inject or task patterns where original form is taken from 
best performer. Comparing 0^3 and 0^3, 0^4 and 0, 4, 
and finally 0^5 with 0^5 shows that small alternation in 
task pattern, obtained by moving one B particles from 
C2 into Ci, lowers performance considerably. Organisms 
with such alternations do not perform well in the network 
like geometry: each of 0^3, o*.4 and o*.s has v > 1. 
Another example, can be obtained from comparison of 
0^3 with o*^- Both inject and task patterns have been 
altered in 0^3. This change is motivated by sequence of 
organisms in first row of table [H] when read from right 
to left. A priori, the organism o».2 could be considered 
to have v < 1, however this is not the case. The actual 
value for v = 418 is (lot) larger than one. 

In summary, the table llTI shows a couple of interesting 
features. First, one can see that when maximum allowed 
number of particles in the system increases new effects 
appear. It is impossible to predict winner in each class 
before calculation is done. Second, there is also a large 
sensitivity on inject and task patterns. Slight alternation 
of these patterns can lead to drastic changes in perfor- 
mance criteria v. 



VII. DISCUSSION 

We introduced what we might call a generic model for 
study of chemical reactions in structured spaces, based 
on a simple way of incorporating interplay between trans- 
port, chemical reactions, and geometry. A number of 
different chemical reactions were studied and their per- 
formance compared in various ways. The main idea is 
to see how the reaction processes behave when geometry 
changes from compact open space towards the cell like 
environment that is more structured. In order to be able 
to analyze large number of reactions, relatively simple 
model was used in order to reduce computational time. 
Following simplifications were made. 

As a model of structured space we used a network 
consisting of containers connected by tubes. Such setup 
captures the most important geometrical aspects of the 
problem. There are compact regions in space that re- 
actants can explore and react within. These regions are 
connected to each other and exchange chemicals. In such 
a way the two most important processes are clearly sepa- 
rated, the local dynamics within container, and the trans- 
port among containers. This is pretty much how living 
cell operates. In such a way one can capture most im- 
portant geometrical aspects of the cell interior. 

The reaction scheme considered is rather crude. Dy- 
namics within container, however complicated it may be, 
is projected onto a relatively low number of degrees of 
freedom, the number of particles in each container. The 
reaction rate A describes how fast reactions happen. The 
description in terms of number of particles becomes more 
and more valid when size of the reaction volume is re- 
duced and number of reactants few 0, 22 • ft is exactly 
this limit that we focusscd on. 

The transport between containers was modeled in sim- 
plest possible way by using the effective transport rate 
D. It was assumed that transport rate is same for any 
pair of containers, and type of particles. This assump- 
tion could be easily relaxed. Clearly, when tubes are 
long, the number of particles in each container decays 
non-exponentially and the transport process can not be 
described by using a transport rate. But such effects are 
not considered here. 

Thus, only two parameters are used to describe dy- 
namics, the reaction rate A and transport rate D. There 
are several advantages in doing so, but the most impor- 
tant one is that reactions and transport are clearly sep- 
arated, and it is easier to understand how they influence 
each other. Despite apparent simplicity, the model stud- 
ied in here captures all essential features of the problem. 
Should there be need for that one can make the model 
more realistic, but we refrain from this at the moment. 

Which types of chemical reactions draw most bene- 
fit from structured spaces? In an attempt to answer 
this question two schemes were formulated. In ROGE 
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scheme one explores variety of chemical reactions while 
spatial structure (geometry) is kept fixed. It is the other 
way around for GORE scheme. At the moment ROGE 
scheme is studied in a lot more detail. The study of 
GORE ensemble will be left for the future work. The 
main findings of this work are discussed bellow. 

Since dynamics is stochastic one needs probabilistic de- 
scription of the system, and in order to perform analysis 
of any chemical reaction one has to solve corresponding 
master equation that describes statistics of events. Two 
techniques for doing this were used. First, the computer 
simulation is a straight forward way to analyze master 
equation (discussed in section ITTTfl , This technique is not 
that accurate due to the stochastic spread of data. One 
needs unrealistically large number of runs in order to gain 
reasonable accuracy in v. In here, a novel method of an- 
alyzing chemical reaction was developed, with emphasis 
on the first passage time (see section Hv|l . This method 
relies on a matrix representation of the master equation. 
In principle, it is exact up to the numerical errors in car- 
rying out matrix manipulations, such as finding inverses, 
multiplying matrix with vector etc. In here we developed 
a equation that describes first passage time distribution 
function for achieving set of tasks. From this equation 
one can easily obtain all moments of distribution function 
and calculate average execution time for reaction and its 
variance. 

It is not so easy to compare two vastly different chem- 
ical reactions. A method for doing this was discussed in 
section[V] There are some obvious difficulties when doing 
comparison. For example, one might need to compare a 
situation where number of molecular types involved in 
reactions are completely different. Also, spatial struc- 
tures need not be the same, e.g. there might be a need 
to compare two vastly different network structures. In 
order to overcome these difficulties one has to use a mea- 
sure of relative performance. For example, in the case of 
a ROGE ensemble the execution time for given reaction 
was compared when geometry is stretched (transport oc- 
curs with rate D n ~ A) and compact (transport with 
rate D where D ^> D n , A, £A). The ratio v of the mag- 
nitude of the execution times t„ and To for stretched 
(network) and compact geometry is a rough estimate of 
how well reaction performs in a network like geometry. 
Instead of comparing organisms directly one compares 
values v = ||r n ||/||To||. 

Thus, the single variable v is sufficient to determine 
whether a reaction (organism) runs best on the network 
like (y < 1) or the compact geometry (y > 1). The 
variable v was used to classify performance of various 
organisms in the geometry consisting of two containers 
connected by tube, with main findings summarized in 
Fig. and Table HU 

(1) It is obvious that intuition does not help much. One 
really has to do numerical analysis in order to extract 
best performer in a given class. For example, the charac- 



ter of best performer changes quite a bit when number of 
particles in the system varies. The structure of the win- 
ning organism in the first row of table changes in a rather 
unpredictable manner as upper limit to the total number 
of particles in the system increases from 2 to 5. Also, the 
performance is extremely sensitive to the details of inject 
and task patterns. It is interesting to speculate whether 
such sensitivity on total particle number, inject and task 
patterns is exploited in the intra cellular machinery. We 
developed a quantitative method to judge on such effects. 

(2) The role of reactions with positive and negative 
catalytic influence is symmetric. Roughly, they occur 
equally often in the organisms that perform well in net- 
work like geometries with v < 1. Interestingly enough, 
same holds for opposite range with v > 1. From Fig. 
one can see that positive and negative areas are roughly 
equal in magnitude, both in the v < 1 and v > 1 re- 
gions. Thus, there are reactions of the suppressor type 
that run better in compact geometries, and the other 
way around. There are also reactions with solely positive 
catalytic influence that run faster in network like geome- 
try. But these findings are the not the only ones that are 
counterintuitive. It is also surprising that in the region 
v ~ reactions with solely positive catalytic influence 
dominate. 

(3) Reaction schemes containing antagonistic catalytic 
influences in the intermediate stages of reaction, that 
slow down the production of the final product, require 
network like geometry to run fast. The antagonistic sub- 
reactions have to isolated and allowed to occur in sep- 
arate regions of space. There are two mechanisms that 
work against each other. First, the reaction time gets 
larger since one needs to transport reactants to the re- 
gions where damage inflicting sub-reactions must be run 
in isolation. Second, after being isolated, dangerous sub- 
reactions happen much faster than when all reactants are 
mixed and this shortens reaction time. It is interesting 
to note that there are large number of cases where isolat- 
ing misbehaving sub-reactions pays off in faster reaction 
time. 

(4) In general, antagonistic catalytic influences are 
hard to identify. It is hard to judge reaction by itself. 
The whole triple consisting of reaction, inject pattern, 
and task pattern has to be considered simultaneously. 

In summary, we studied the workings of a chemical 
reactions in a two vastly different spaces. In the com- 
pact conformation, molecules can reach any part of the 
space very fast. The transport rate is lot larger than 
any reaction rate in the system. In the network like con- 
formation, various volumes are well separated and the 
transport of reactants between volumes occurs relatively 
slowly. It was found that considerable number of reac- 
tions work better in a network like configuration when 
transport rate gets smaller. The setup suggested in here 
is generic. There are many possible ways of extending 
present analysis. For example, at present focus is on 
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small number of particles and stochastic dynamics. One 
could easily consider situation where number of particles 
is larger and classical chemical kinetics applies in the con- 
tainer. Transport between containers can be treated in a 
better way. Instead of focusing on the average execution 
time r one can easily look at noise er or combination of 
the two. The GORE ensemble should be explored in a lot 
more detail. These issues will be addressed in the future 
work. 
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Figures 



(a) 



(1)) 



(c) 




FIG. 1: A guidance how to think about the problem. Panel 

(a) : non-structured (compact) space. Panel (b): structured 
space obtained by deforming geometry depicted in panel (a). 
The goal is to understand how reaction dynamics alters when 
geometry changes from (a) to (b). The panel (c), represent- 
ing a network of containers and tubes, is used to achieve this 
goal. For example, to capture the most important geometri- 
cal features of the structure in panel (b) one needs five con- 
tainers Ci,...,Cg and four tubes with lengths Z15, £25, I35 
and I45. The two parameters are used to define the dynam- 
ics in the network (c). The intra-container reaction rate A, 
and inter-container transport rate D (provided containers are 
connected). When D ~ A the reaction dynamics in panels 

(b) and (c) should exhibit some similarities. For D>) one 
expects the same for the structures (a) and (c). 




FIG. 2: The most general form of a reaction graph consid- 
ered. Species Xi,X<z, . . . , X u describe an inject pattern, a set 
of the molecules that are inserted into the system at t = 0. 
This set is arranged into the vector t. Molecules denoted by 
X[ , X'2, • • • , X' v are the ones that have to be synthesized and 
are arranged into the vector 7r . These are referred to as a task 
pattern. The shaded area in the middle denotes intermediate 
reaction steps. Dashed lines with arrows are allowed and re- 
sult in appearance of loops in the reaction scheme. The speed 
of reaction is calculated by tracking instances when particles 
in the set tt appear for the first time. The corresponding times 
T\, T%, . . . , T v are arranged into the vector T. See section [D] 
for more details. 
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§ t - trap states ; 



FIG. 3: The structure of the configuration space. Initially the 
system is in the state i. Three set of states are distinguished, 
a set of normal states S n , a set of trap states St, and a set of 
window states S m . See text for discussion. 
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FIG. 5: Classification of organisms in ROGE ensemble, v was 
calculated with D n — Is -1 , Do = 3125s -1 , A = Is" 1 , and 
£ = 100. Panels (a), (c), (e), (g) and (i) are histograms that 
depicts groups with similar reaction speed v. There are 100 
classes of width 0.01 for v from 1 to 2. For v > 1, instead 
of v, the value obtained from x( u ) = 2 — ln2/ln(l + v) is 
used. The function \ monotonically grows with v and maps 
infinite interval [1, oo) onto the finite one [1,2). In addition, 
this particular form for x reveals more details in the region 
near v = 1. Panels (b), (d), (f), (h), and (j) are discrete 
spectra (no histogram) for region v £ [0,1]. Panels in the 
same row have same value for the total number of particles 
in the system N*: (a) and (b) N* = 1, (c) and (d) TV* = 2, 
(e) and (f) N; = 3, (g) and (h) iV p * = 4, (i) and (j) 7V p * = 
5. Negative value for number of organisms (#o) indicates 
that organism with particular value of v contains at least 
one reaction that can be influenced through mechanism of 
suppression. Presence of suppressor lowers the reaction rate 
from A to A/5 (see sect ion HT1 for details). 



FIG. 4: The dependence of r and a on D -1 ~ yfl for three 
organisms defined in tableQ](Z is the length of a tube connect- 
ing containers). Curves describe oi (dotted line), 02 (dashed) 
and 03 (solid). See table and section IVI for more details. 
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TABLE I: Analysis of a three cases where the total number 
of particles equals one and the number of containers equals 
two. These are simplest cases that one can consider with 
regard to the choices of l and n. A reaction scheme used is 
A — > B and B — + A, both with rate A (effects of catalyst or 
suppressor are absent since there is only one particle in the 
system at a time). The jump rate between containers is given 
by D. In the definition of i and tt the content of containers 
is symbolically indicated as [content of Ci]-[content of C2]. 



org. 
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TT 
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a 2 


01 
02 
03 


[A]-[ ] 
[A]-[ ] 
[A]-[ ] 


[ ]-[A] 
[B]-[ ] 
[ ]-[B] 


2 D+X 


2 ZJ 3 +5 D 2 X+3D X 2 + X s 


D(D+X) 
D+2 X 


D 2 X (D + X)' 2 
D 3 +3 D 2 A+5 D X 2 + 2 A 3 


X (D+X) 

- + - 

D X 


D A 2 (D + X) 2 
DX 



TABLE II: List of the best performers in the ROGE ensemble (first row). Couple of second best performers are also shown 
(second row) . Last row contains organisms that are obtained by slightly perturbing winners from the first row. 



n; = 2 


n; = 3 


n; = 4 


n; = 5 


( Oi,2 ) 

A —-% B, B A (iii) 


( °i,3 ) 

A M B (R 3 ) 


( Oi,4 ) 

a±Ab (R 3 ) 


( Oi,8 ) 

A ^ B (i? 3 ) 


t =[A]-[B] 


t =[A]-[2A] 


l =[A]-[3A] 


i =[2A]-[3A] 


7T =[2B]-[ ] or [ ]-[2B] 


7r=[2A]-[B] 


tt =[3A]-[B] 


tt =[4A]-[B] 


t„=[18.55]-[ ] 


t„=[0.242]-[0.0355] 


t„ =[0.021]-[0.01] 


r n =[0.014]-[0.00063] 


r =[50.91]-[ ] 


t =[0.000219]-[0.979] 


r =[0.00045]-[0.95] 


r =[0.010]-[0.91] 


v = 0.364 


^=0.250 


= 0.0248 


v = 0.0192 


( °«,2 ) 


( °«,3 ) 






A — % B, B -±4 A (R 2 ) 


Jfc, i=[ ]-[3A] 


( Oii,4 ) 

i? 3 , i=[2A]-[2A] 
tt=[3A]-[B], 1/ =0.0279; 

Rs, t=[ ]-[4A] 
tt=[3A]-[B], i/=0.0358 


( 0«,5 ) 

Rs, t=[2A]-[3A] 
tt=[4A,B]-[ ], ^=0.0219; 
Rs, t=[A]-[4A] 
tt=[4A]-[B], !/=0.0305 


1 =[A]-[B], tt =[2B]-[ ] 
or [ ]-[2B], v = 0.368 

Ri,l=[ ]-[2A], 
7T = [2B]-[ ],v = 0.386 

Rs,i=[ ]-[2A], 


tt=[2A]-[B], i/ = 0.268 ; 

Ri, t=[2A]-[B] 
tt=[A]-[2B] 

R 2 , t=[2A]-[B] 
w=[2B]-[A] or [ ]-[A,2B] 


7T =[A]-[B], 1/ = 0.500 


v « 0.386 






( °*,2 ) 


( °*,3 ) 


( °*A ) 


( °*,5 ) 


i?i or R 3 , 1 =[ ]-[2A] 


Rs, 1 =[A]-[2A] 


Jfe, i =[A]-[3A] 


Rs, 1 =[2A]-[3A] 


7T =[2A]-[ ], i/ =3249 or 418; 


7T =[3A]-[ \,v = 21.8 


tt =[4A]-[ ], 1/= 17.5 


tt =[5A]-[ ],v= 10.4 



